function comoving_distance_integrand(z) result(integrand)
  real*4 :: z
  real*4 :: integrand
  real*4 :: omega0_common, lambda0_common
  common /comoving_common/ omega0_common, lambda0_common

  integrand = 1.0e0/sqrt((1.0e0+z)**2.0e0*(1.0e0+omega0_common*z)-z*(2.0e0+z)*lambda0_common)
  return

end function comoving_distance_integrand

function compute_comoving_distance(z,omega0,lambda0,hubble0) result(comdist)
  implicit none
  real*4 :: integral
  real*4 :: omega0_common, lambda0_common
  real*8 :: z,omega0,lambda0,hubble0
  real*8 :: omegak
  real*8 :: comdist
  common /comoving_common/ omega0_common, lambda0_common

  if(lambda0 < 1.0d-10)then

     comdist = omega0*z + (omega0-2.0d0)*(dsqrt(1.0d0+omega0*z)-1.0d0)
     comdist = comdist/omega0**2.0d0/(1.0d0+z)
     comdist = comdist * 2.0d0 * 2.9979d3/hubble0

  else

     ! ** Carroll & Press 1992, Eqn. 25 **
     omega0_common  = real(omega0)
     lambda0_common = real(lambda0)
     call my_qsimp('comoving_distance_integrand',0.0e0,real(z),integral)
   
     omegak = 1.0d0 - omega0 - lambda0

     if(abs(omegak) < 1.0e-16) then 
        comdist = dble(integral)
     else
        if(omegak > 0.0e0)then
           comdist = dsinh(dsqrt(abs(omegak))*dble(integral))
        else 
           comdist = dsin(dsqrt(abs(omegak))*dble(integral))
        endif
        comdist = comdist / dsqrt(abs(omegak))
     endif
     comdist = comdist * 2.9979d3/hubble0

  endif

end function compute_comoving_distance
